% By Martin M. Andreasen, Nov 2016
% This function simulates the model when solved up to third order. 
% The purning scheme is NOT used and we apply the kronecker notation.
% The simulated state space system reads:
% 
% y_t   = g(x_t,sig)
% x_t+1 = h(x_t,sig) + sig*eta*shocks_t+1 
%
function [Y_sim,X_sim] = sim3rdNoPrun(model,shocks,orderApp)

% Unfold model
gx  = model.gx;
gxx = model.gxx;
gxxx= model.gxxx;
gss = model.gss;
gssx= model.gssx;
gsss= model.gsss;

hx  = model.hx;
hxx = model.hxx;
hxxx= model.hxxx;
hss = model.hss;
hssx= model.hssx;
hsss= model.hsss;

eta = model.eta;
sig = 1;

numSim = size(shocks,2);    %Number of simulations

% Allocating memory
Y_sim   = zeros(model.ny,numSim);
X_sim   = zeros(model.nx,numSim);

% The simulation
xt    = zeros(model.nx,1);    %The first order terms

% Defining matrices
HHxxtil    = 1/2*reshape(hxx,model.nx,model.nx^2);
HHxxxtil   = 1/6*reshape(hxxx,model.nx,model.nx^3);
GGxxtil    = 1/2*reshape(gxx,model.ny,model.nx^2);   
GGxxxtil   = 1/6*reshape(gxxx,model.ny,model.nx^3);   
X_sim(:,1) = xt;
if orderApp == 1
   for t=1:numSim
      X_sim(:,t) = xt; 
      Y_sim(:,t) = gx*xt;
       
      % The states next period
      xt_p  = hx*xt + sig*eta*shocks(:,t);

      % Updating xt
      xt  = xt_p;
   end
elseif orderApp == 2
   for t=1:numSim
      X_sim(:,t) = xt;
      AA         = kron(xt,xt);
      Y_sim(:,t) = gx*xt + GGxxtil*AA + 1/2*sig^2*gss; 
       
      % The states mext period
      xt_p  = hx*xt + sig*eta*shocks(:,t) + HHxxtil*AA + 1/2*sig^2*hss;
    
      % Updating xt
      xt  = xt_p;
   end
elseif orderApp == 3
   for t=1:numSim
      X_sim(:,t) = xt;
      AA    = kron(xt,xt);
      AAA   = kron(xt,AA);
      Y_sim(:,t) = gx*xt + GGxxtil*AA + 1/2*sig^2*gss + GGxxxtil*AAA + 3/6*gssx*sig^2*xt + 1/6*sig^3*gsss;
      
      % The states next period
      xt_p  = hx*xt + sig*eta*shocks(:,t) + HHxxtil*AA + 1/2*sig^2*hss + HHxxxtil*AAA+ 3/6*hssx*sig^2*xt + 1/6*sig^3*hsss;
      
      % Updating xt
      xt = xt_p;
   end
end
